setwd("D:/Valdivia/DissertationChapter3/RCode/")


library(ggpmisc)
library(EcoHydRology)
library(tidyverse)
library(tidyr)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(sf)
library(gridExtra)
library(broom)

#Creating table to show how flooding at nodes have changed between the base model and the modified model in the subsystem
comparisonData <- data.frame(read.csv("BaseModifiedComparison_PresDay.csv"))
comparisonData <- comparisonData[which(comparisonData$FloodVolume > 0), ]

uniqueNodes_Base <- unique(comparisonData$Node[which(comparisonData$Model == "Base")])
uniqueNodes_Modified <- unique(comparisonData$Node[which(comparisonData$Model == "Modified")])

comparisonData_improved <- data.frame(matrix(nrow = 0, ncol = 3, NA))
colnames(comparisonData_improved) <- c("Node", "FloodVolume_Base", "FloodVolume_Modified")
i = 1
for(i in 1:length(uniqueNodes_Base)){
  rowToAdd <- data.frame(matrix(nrow = 1, ncol = 3, NA))
  colnames(rowToAdd) <- c("Node", "FloodVolume_Base", "FloodVolume_Modified")
  
  rowToAdd[1, 1] <- uniqueNodes_Base[i]
  rowToAdd[1, 2] <- comparisonData$FloodVolume[which(comparisonData$Node == uniqueNodes_Base[i] & comparisonData$Model == "Base")]
  if(length(comparisonData$FloodVolume[which(comparisonData$Node == uniqueNodes_Base[i] & comparisonData$Model == "Modified")]) > 0){
    rowToAdd[1, 3] <- comparisonData$FloodVolume[which(comparisonData$Node == uniqueNodes_Base[i] & comparisonData$Model == "Modified")]
  }
  comparisonData_improved <- rbind(comparisonData_improved, rowToAdd)
}

for(i in 1:length(uniqueNodes_Modified)){
  if(length(which(uniqueNodes_Base == uniqueNodes_Modified[i])) == 0){
    rowToAdd <- data.frame(matrix(nrow = 1, ncol = 3, NA))
    colnames(rowToAdd) <- c("Node", "FloodVolume_Base", "FloodVolume_Modified")
    
    rowToAdd[1, 1] <- uniqueNodes_Modified[i]
    rowToAdd[1, 3] <- comparisonData$FloodVolume[which(comparisonData$Node == uniqueNodes_Modified[i] & comparisonData$Model == "Modified")]
    
    comparisonData_improved <- rbind(comparisonData_improved, rowToAdd)
  }
}

comparisonData_improved$FloodVolume_Base[is.na(comparisonData_improved$FloodVolume_Base)] <- 0
comparisonData_improved$FloodVolume_Modified[is.na(comparisonData_improved$FloodVolume_Modified)] <- 0
comparisonData_improved$Difference_FloodVolume <- comparisonData_improved$FloodVolume_Modified - comparisonData_improved$FloodVolume_Base

length(which(comparisonData_improved$Difference_FloodVolume < -0.01))
length(which(comparisonData_improved$Difference_FloodVolume > 0.01))
length(which(comparisonData_improved$Difference_FloodVolume <= 0.01 & comparisonData_improved$Difference_FloodVolume >= -0.01))
sum(abs(comparisonData_improved$Difference_FloodVolume))

sum(comparisonData_improved$FloodVolume_Base[!is.na(comparisonData_improved$FloodVolume_Base)])
sum(comparisonData_improved$FloodVolume_Modified[!is.na(comparisonData_improved$FloodVolume_Modified)])

length(comparisonData_improved$FloodVolume_Base[!is.na(comparisonData_improved$FloodVolume_Base)])
length(comparisonData_improved$FloodVolume_Modified[!is.na(comparisonData_improved$FloodVolume_Modified)])


#Downstream node flooding analysis
wet43_maxDepth <- 1.25
S11_01_maxDepth <- 1.61
S11_05_maxDepth <- 0.88
S11_06_maxDepth <- 1.45
scenarios <- c("BAU", "EcoWet", "Friendly", "Inclusive", "R2F")
modelTypes <- c("Wetland", "Wshed", "WetlandWshed")
nodes <- c("wet43", "S11_01", "S11_05", "S11_06")
conduits <- c("ConduitWe43", "CuS11_01", "CuS11_04", "EmS11_05", "EmS11_06", "CuS11_05", "CuS11_06a", "CuS11_06b")

#Use these for tabulating the Present Day scenario output info
#scenarios <- c("PresDay")
modelType <- c("Wetland")


i = 1
j = 1
k = 1


#This script tabulates the flooding at each node in the Wetland43 system of the EPA Swmm. The only required inputs are the .txt table outputs of the node depths and the conduit flow data. After exporting a .txt from EPA SWMM, you need to go and delete the header off of all the text files manually.
for(i in 1:length(scenarios)){
  for(j in 1:length(modelType)){
    
    eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_allNodes <- data.frame(matrix(nrow = 0, ncol = 6, NA))", sep = "")))
    eval(parse(text = paste("colnames(floodVolumes_", scenarios[i], "_", modelType[j], "_allNodes) <- c(\"Scenario\", \"Strength\", \"WetlandLoss\", \"NodeFlooding\", \"Feature\", \"Node\")", sep = "")))
    
    
    for(k in 1:length(nodes)){
      eval(parse(text = paste(nodes[k], "_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/", nodes[k], "_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
      eval(parse(text = paste("colnames(", nodes[k], "_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Depth\")", sep = "")))
      eval(parse(text = paste(nodes[k], "_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
      
      eval(parse(text = paste("floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], " <- ", nodes[k], "_", scenarios[i], "_", modelType[j], "[which(", nodes[k], "_", scenarios[i], "_", modelType[j], "$Depth >= ", nodes[k], "_maxDepth), 1:2]", sep = "")))
    }
    
    for(k in 1:length(nodes)){
    
      if(k == 1) {
        eval(parse(text = paste("conduitWet43_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/ConduitWet43_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(conduitWet43_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("conduitWet43_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(conduitWet43_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("conduitWet43_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(conduitWet43_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        
        floodVolumes_wet43 <- data.frame(matrix(nrow = 0, ncol = 3, NA))
        colnames(floodVolumes_wet43) <- c("Days", "Time", "Volume_m3")
        
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- conduitWet43_", scenarios[i], "_", modelType[j], "[which(conduitWet43_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & conduitWet43_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_wet43 <- rbind(floodVolumes_wet43, floodRow)}}\n}", sep = "")))
        
        floodVolumes_wet43$Flow_cms <- floodVolumes_wet43$Flow_cms * 60 / 1000
        floodVolumes_wet43$Flow_cms[which(floodVolumes_wet43$Flow_cms < 0)] <- floodVolumes_wet43$Flow_cms[which(floodVolumes_wet43$Flow_cms < 0)] * -1
        }
      
          
    
      if(k == 2){
        eval(parse(text = paste("CuS11_01_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/CuS11_01_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(CuS11_01_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("CuS11_01_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(CuS11_01_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("CuS11_01_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(CuS11_01_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        
        floodVolumes_S11_01 <- data.frame(matrix(nrow = 0, ncol = 3, NA))
        colnames(floodVolumes_S11_01) <- c("Days", "Time", "Volume_m3")
        
        
        
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- conduitWet43_", scenarios[i], "_", modelType[j], "[which(conduitWet43_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & conduitWet43_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms < 0){floodVolumes_S11_01 <- rbind(floodVolumes_S11_01, floodRow)}}\n}", sep = "")))
        
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- CuS11_01_", scenarios[i], "_", modelType[j], "[which(CuS11_01_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & CuS11_01_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_S11_01 <- rbind(floodVolumes_S11_01, floodRow)}}\n}", sep = "")))
        
        
        floodVolumes_S11_01$Flow_cms <- floodVolumes_S11_01$Flow_cms * 60 / 1000
        floodVolumes_S11_01$Flow_cms[which(floodVolumes_S11_01$Flow_cms < 0)] <- floodVolumes_S11_01$Flow_cms[which(floodVolumes_S11_01$Flow_cms < 0)] * -1
        
      } 
      

      if(k == 3){
        eval(parse(text = paste("CuS11_04_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/CuS11_04_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(CuS11_04_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("CuS11_04_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(CuS11_04_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("CuS11_04_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(CuS11_04_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        eval(parse(text = paste("EmS11_05_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/EmS11_05_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(EmS11_05_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("EmS11_05_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(EmS11_05_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("EmS11_05_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(EmS11_05_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        eval(parse(text = paste("CuS11_05_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/CuS11_05_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(CuS11_05_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("CuS11_05_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(CuS11_05_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("CuS11_05_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(CuS11_05_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        
        
        floodVolumes_S11_05 <- data.frame(matrix(nrow = 0, ncol = 3, NA))
        colnames(floodVolumes_S11_05) <- c("Days", "Time", "Volume_m3")
        
        
        #For updrain conduit CuS11_04. Flow should be negative for it to be from S11_05 because inlet node is S11_04 and outlet node is S11_05
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- CuS11_04_", scenarios[i], "_", modelType[j], "[which(CuS11_04_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & CuS11_04_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms < 0){floodVolumes_S11_05 <- rbind(floodVolumes_S11_05, floodRow)}}\n}", sep = "")))
        
        #For embalse EMB16. Flow should be positive because inlet node is S11_05 and outlet is EMB16
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- EmS11_05_", scenarios[i], "_", modelType[j], "[which(EmS11_05_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & EmS11_05_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_S11_05 <- rbind(floodVolumes_S11_05, floodRow)}}\n}", sep = "")))
        
        #For downdrain conduit CuS11_05. Flow should be positive
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- CuS11_05_", scenarios[i], "_", modelType[j], "[which(CuS11_05_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & CuS11_05_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_S11_05 <- rbind(floodVolumes_S11_05, floodRow)}}\n}", sep = "")))
        
        floodVolumes_S11_05$Flow_cms <- floodVolumes_S11_05$Flow_cms * 60 / 1000
        floodVolumes_S11_05$Flow_cms[which(floodVolumes_S11_05$Flow_cms < 0)] <- floodVolumes_S11_05$Flow_cms[which(floodVolumes_S11_05$Flow_cms < 0)] * -1
        
      }
      
      if(k == 4){
        eval(parse(text = paste("EmS11_06_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/EmS11_06_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(EmS11_06_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("EmS11_06_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(EmS11_06_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("EmS11_06_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(EmS11_06_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        eval(parse(text = paste("CuS11_06a_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/CuS11_06a_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(CuS11_06a_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("CuS11_06a_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(CuS11_06a_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("CuS11_06a_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(CuS11_06a_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        eval(parse(text = paste("CuS11_06b_", scenarios[i], "_", modelType[j], " <- data.frame(read.table(\"D:/Valdivia/DissertationChapter3/RCode/DowndrainOutputs/CuS11_06b_", scenarios[i], "_", modelType[j], ".txt\", skip = 4))", sep = "")))
        eval(parse(text = paste("colnames(CuS11_06b_", scenarios[i], "_", modelType[j], ") <- c(\"Days\", \"Time\", \"Flow_cms\")", sep = "")))
        eval(parse(text = paste("CuS11_06b_", scenarios[i], "_", modelType[j], "$Time <- as.POSIXct(CuS11_06b_", scenarios[i], "_", modelType[j], "$Time, format = \"%H:%M:%S\")", sep = "")))
        eval(parse(text = paste("CuS11_06b_", scenarios[i], "_", modelType[j], "$Flow_cms <- as.numeric(CuS11_06b_", scenarios[i], "_", modelType[j], "$Flow_cms)", sep = "")))
        
        
        
        floodVolumes_S11_06 <- data.frame(matrix(nrow = 0, ncol = 3, NA))
        colnames(floodVolumes_S11_06) <- c("Days", "Time", "Volume_m3")
        
        #For updrain conduit CuS11_05. Flow should be negative for it to be from S11_06 because inlet node is S11_05 and outlet node is S11_06
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- CuS11_05_", scenarios[i], "_", modelType[j], "[which(CuS11_05_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & CuS11_05_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms < 0){floodVolumes_S11_06 <- rbind(floodVolumes_S11_06, floodRow)}}\n}", sep = "")))
        
        #For embalse EMB16. Flow should be positive because inlet node is S11_06 and outlet is EMB16
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- EmS11_06_", scenarios[i], "_", modelType[j], "[which(EmS11_06_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & EmS11_06_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_S11_06 <- rbind(floodVolumes_S11_06, floodRow)}}\n}", sep = "")))
        
        #For downdrain conduit CuS11_06a. Flow should be positive
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- CuS11_06a_", scenarios[i], "_", modelType[j], "[which(CuS11_06a_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & CuS11_06a_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_S11_06 <- rbind(floodVolumes_S11_06, floodRow)}}\n}", sep = "")))
        
        #For downdrain conduit CuS11_06b. Flow should be positive
        eval(parse(text = paste("if(length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time) > 0){\nfor(l in 1:length(floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time)){\nfloodRow <- CuS11_06b_", scenarios[i], "_", modelType[j], "[which(CuS11_06b_", scenarios[i], "_", modelType[j], "$Days == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Days[l] & CuS11_06b_", scenarios[i], "_", modelType[j], "$Time == floodTimes_", nodes[k], "_", scenarios[i], "_", modelType[j], "$Time[l]),] \nif (floodRow$Flow_cms > 0){floodVolumes_S11_06 <- rbind(floodVolumes_S11_06, floodRow)}}\n}", sep = "")))
        
        
        floodVolumes_S11_06$Flow_cms <- floodVolumes_S11_06$Flow_cms * 60 / 1000
        floodVolumes_S11_06$Flow_cms[which(floodVolumes_S11_06$Flow_cms < 0)] <- floodVolumes_S11_06$Flow_cms[which(floodVolumes_S11_06$Flow_cms < 0)] * -1
        
      }
      
      #Tabulating results for all nodes for a given scenario and model type
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], " <- data.frame(matrix(nrow = 1, ncol = 6, NA))", sep = "")))
      eval(parse(text = paste("colnames(floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], ") <- c(\"Scenario\", \"Strength\", \"WetlandLoss\", \"NodeFlooding\", \"Feature\", \"Node\")", sep = "")))
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], "[1,1] <- \"", scenarios[i], "\"", sep = "")))
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], "[1,2] <- NA", sep = "")))
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], "[1,3] <- 0", sep = "")))
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], "[1,4] <- sum(floodVolumes_", nodes[k], "$Flow_cms)", sep = "")))
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], "[1,5] <- \"", modelType[j], "\"", sep = "")))
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], "[1,6] <- \"", nodes[k], "\"", sep = "")))
      
      eval(parse(text = paste("floodVolumes_", scenarios[i], "_", modelType[j], "_allNodes <- rbind(floodVolumes_", scenarios[i], "_", modelType[j], "_allNodes, floodVolumes_", scenarios[i], "_", modelType[j], "_", nodes[k], ")", sep = "")))
      
    }
    
    

  }
}




scenarios <- c("BAU", "EcoWet", "Friendly", "Inclusive", "R2F")
modelType <- c("Wetland", "Wshed", "WetlandWshed")
nodes <- c("wet43", "S11_01", "S11_05", "S11_06")
conduits <- c("ConduitWe43", "CuS11_01", "CuS11_04", "EmS11_05", "EmS11_06", "CuS11_05", "CuS11_06a", "CuS11_06b")


eval(parse(text = paste("floodVolumes_Master <- data.frame(matrix(nrow = 0, ncol = 6, NA))", sep = "")))
eval(parse(text = paste("colnames(floodVolumes_Master) <- c(\"Scenario\", \"Strength\", \"WetlandLoss\", \"NodeFlooding\", \"Feature\", \"Node\")", sep = "")))
  
#Use these for tabulating the Present Day scenario output info
scenarios <- c("PresDay")
modelType <- c("Wetland")



for(i in 1:length(scenarios)){
  eval(parse(text = paste("floodVolumes_Master <- rbind(floodVolumes_Master, floodVolumes_", scenarios[i], "_", modelType, "_allNodes)", sep = "")))
}


for(i in 1:length(floodVolumes_Master$Scenario)){
  if(floodVolumes_Master$Scenario[i] == "BAU"){
    floodVolumes_Master$Scenario[i] <- "Business-as-usual"
    floodVolumes_Master$Strength[i] <- "Low"
    floodVolumes_Master$WetlandLoss[i] <- 0.97
  } else if(floodVolumes_Master$Scenario[i] == "EcoWet"){
    floodVolumes_Master$Scenario[i] <- "Eco-wetland"
    floodVolumes_Master$Strength[i] <- "High"
    floodVolumes_Master$WetlandLoss[i] <- 0.76
  } else if (floodVolumes_Master$Scenario[i] == "Friendly") {
    floodVolumes_Master$Strength[i] <- "Med"
    floodVolumes_Master$WetlandLoss[i] <- 0.30
  } else if (floodVolumes_Master$Scenario[i] == "Inclusive") {
    floodVolumes_Master$Strength[i] <- "Med"
    floodVolumes_Master$WetlandLoss[i] <- 0.35
  } else if(floodVolumes_Master$Scenario[i] == "R2F"){
    floodVolumes_Master$Scenario[i] <- "Resilient to flooding"
    floodVolumes_Master$Strength[i] <- "High"
    floodVolumes_Master$WetlandLoss[i] <- 0.66
  } else if(floodVolumes_Master$Scenario[i] == "PresDay"){
    floodVolumes_Master$Scenario[i] <- "Present day"
    floodVolumes_Master$WetlandLoss[i] <- 0.00
  }
}

presDay_WetlandWshed <- floodVolumes_Master[which(floodVolumes_Master$Scenario == "Present day"),]
presDay_WetlandWshed$Feature <- "WetlandWshed"
presDay_Wshed <- floodVolumes_Master[which(floodVolumes_Master$Scenario == "Present day"),]
presDay_Wshed$Feature <- "Wshed"

floodVolumes_Master <- rbind(floodVolumes_Master, presDay_WetlandWshed, presDay_Wshed)

#write.csv(floodVolumes_Master, file = "floodVolumes_Master_Wetland43andDowndrain.csv")

#Reading the file created in the previous step, as a shortcut to doing all the calculations again
floodVolumes_Master <- data.frame(read.csv("floodVolumes_Master_Wetland43andDowndrain.csv"))

i = 1
j = 1
k = 1
#Collecting significance for all downdrain flooding nodes in subsystem 1
##Building variables and data frame

master_pValues_df <- data.frame(matrix(ncol = 4, nrow = 0, NA))
colnames(master_pValues_df) <- c("Node", "modelType", "pValue", "rSquared")
for(k in length(nodes)){

  eval(parse(text = paste(nodes[k], "_pValues_df <- data.frame(matrix(ncol = 4, nrow = 3, NA))", sep = "")))
  eval(parse(text = paste("colnames(", nodes[k], "_pValues_df) <- c(\"Node\", \"modelType\", \"pValue\", \"rSquared\")", sep = "")))
  eval(parse(text = paste(nodes[k], "_pValues_df$Node <- \"", nodes[k], "\"", sep = "")))
  eval(parse(text = paste("floodVolumes_", nodes[k], " <- floodVolumes_Master[which(floodVolumes_Master$Node == \"", nodes[k], "\"),]", sep = "")))
  
  for(j in 1:length(modelType)){
    eval(parse(text = paste(nodes[k], "_pValues_df[j, 2] <- \"", modelType[j], "\"", sep = "")))
    eval(parse(text = paste("lm_", nodes[k], "_", modelType[j], " <- lm(NodeFlooding ~ WetlandLoss, floodVolumes_", nodes[k], "[which(floodVolumes_", nodes[k], "$Feature == \"", modelType[j], "\"),])", sep = "")))
    eval(parse(text = paste(nodes[k], "_pValues_df[j, 3] <- glance(lm_", nodes[k], "_", modelType[j], ")$p.value", sep = "")))
    eval(parse(text = paste(nodes[k], "_pValues_df[j, 4] <- summary(lm_", nodes[k], "_", modelType[j], ")$adj.r.squared", sep = "")))
  }
}
eval(parse(text = paste("master_pValues_df <- rbind(master_pValues_df, ", nodes, "_pValues_df)", sep = "")))


#write.csv(master_pValues_df, file = "master_pValues_df_Wet43AndDowndrain.csv")
master_pValues_df <- read.csv("master_pValues_df_Wet43AndDowndrain.csv")



#Plotting flooding at downdrain nodes in subsystem 1
eval(parse(text = paste("g_Wet43 <- ggplot(data = floodVolumes_Master[which(floodVolumes_Master$Node == \"wet43\"),], aes(x = WetlandLoss, y = NodeFlooding, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))

eval(parse(text = paste("g_Wet43_1 <- g_Wet43 + labs(title = \"Updrain wetland\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"solid\", \"blank\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1)) +
                        scale_y_continuous(limits = c(0, 12), breaks = seq(0, 12, 2))",
                        sep = "")))



eval(parse(text = paste("g_S11_01 <- ggplot(data = floodVolumes_Master[which(floodVolumes_Master$Node == \"S11_01\"),], aes(x = WetlandLoss, y = NodeFlooding, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))


eval(parse(text = paste("g_S11_01_1 <- g_S11_01 + labs(title = \"Downdrain node 1\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost updrain wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"solid\", \"dotted\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1)) + 
                        scale_y_continuous(limits = c(0, 12), breaks = seq(0, 12, 2))",
                        sep = "")))



eval(parse(text = paste("g_S11_05 <- ggplot(data = floodVolumes_Master[which(floodVolumes_Master$Node == \"S11_05\"),], aes(x = WetlandLoss, y = NodeFlooding, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))


eval(parse(text = paste("g_S11_05_1 <- g_S11_05 + labs(title = \"Downdrain node 2\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost updrain wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"solid\", \"blank\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1)) +
                        scale_y_continuous(limits = c(0, 12), breaks = seq(0, 12, 2))",
                        sep = "")))



eval(parse(text = paste("g_S11_06 <- ggplot(data = floodVolumes_Master[which(floodVolumes_Master$Node == \"S11_06\"),], aes(x = WetlandLoss, y = NodeFlooding, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))


eval(parse(text = paste("g_S11_06_1 <- g_S11_06 + labs(title = \"Downdrain node 3\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost updrain wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"blank\", \"blank\", \"blank\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1)) +
                        scale_y_continuous(limits = c(0, 12), breaks = seq(0, 12, 2))",
                        sep = "")))


#Creating a pdf of downdrain flooding at subsystem 1
pdf(file = "Wetland43_LinearPlot_09152021.pdf")
ggarrange(g_Wet43_1, g_Wet43_1,  
          g_S11_05_1, g_S11_06_1,
          ncol = 2, nrow =2)
ggarrange(g_S11_01_1, g_S11_01_1, g_S11_01_1, g_S11_01_1,
          ncol = 2, nrow = 2)
dev.off()


i = 1
j = 1
k = 1

#Updating flooded node shapefile to include node flood volumes and durations for all scenarios and present day
##Scenario names and model types
scenarioNames <- c("Friendly", "Inclusive", "R2F", "EcoWet", "BAU")
modelTypes <- c("Wetlands", "Wsheds", "Wetlands_Wsheds")

##Preparing dataframe of node flooding
for(i in 1:length(modelTypes)){
  for(j in 1:length(scenarioNames)){
    eval(parse(text = paste("nodeFlooding_", modelTypes[i], "_", scenarioNames[j], " <- data.frame(read.csv(\"", modelTypes[i], "_", scenarioNames[j], ".csv\"))", sep = "")))
    eval(parse(text = paste("nodeFlooding_", modelTypes[i], "_", scenarioNames[j], " <- nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$VolumeFlooded_1000m3 >= 0.001),]", sep = "")))
    eval(parse(text = paste("colnames(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], ") <- c(\"Node\", \"FloodDuration_", modelTypes[i], "_", scenarioNames[j], "\", \"FloodVolume_", modelTypes[i], "_", scenarioNames[j], "\")", sep = "")))
  }
}

##Reading table of flooded node locations from shapefile. Note: these locations in the shapefile were created by hand using node flooding information from the SWMM reports and relating them to locations in ArcGIS (also by hand)
floodedNodesMaster <- st_read("FloodedNodes.shp")
floodedNodesMaster <- floodedNodesMaster[,1:2]
nodeList <- floodedNodesMaster$Node

##Adding flood volume and duration information based on secnario and model
for(i in 1:length(modelTypes)){
  for(j in 1:length(scenarioNames)){
    eval(parse(text = paste("floodedNodesMaster$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "<- 0", sep = "")))
    eval(parse(text = paste("floodedNodesMaster$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "<- 0", sep = "")))
    
    for(k in 1:length(nodeList)){
      eval(parse(text = paste("if(length(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]) > 0){floodedNodesMaster$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "[which(floodedNodesMaster$Node == \"", nodeList[k], "\")] <- nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]}", sep = "")))
      eval(parse(text = paste("if(length(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]) > 0){floodedNodesMaster$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "[which(floodedNodesMaster$Node == \"", nodeList[k], "\")] <- nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]}", sep = "")))
    }
  }
}



##PresDay must be done separately because there is only one model for it
floodedNodesMaster$FloodDuration_Wetlands_Wsheds_PresDay<- 0
floodedNodesMaster$FloodVolume_Wetlands_Wsheds_PresDay<- 0

modelTypes <- "Wetlands_Wsheds"
scenarioNames <- "PresDay"
i = 1
j = 1
k = 1

for(i in 1:length(modelTypes)){
  for(j in 1:length(scenarioNames)){
    eval(parse(text = paste("nodeFlooding_", modelTypes[i], "_", scenarioNames[j], " <- data.frame(read.csv(\"", modelTypes[i], "_", scenarioNames[j], ".csv\"))", sep = "")))
    eval(parse(text = paste("nodeFlooding_", modelTypes[i], "_", scenarioNames[j], " <- nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$VolumeFlooded_1000m3 >= 0.001),]", sep = "")))
    eval(parse(text = paste("colnames(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], ") <- c(\"Node\", \"FloodDuration_", modelTypes[i], "_", scenarioNames[j], "\", \"FloodVolume_", modelTypes[i], "_", scenarioNames[j], "\")", sep = "")))
  }
}

for(k in 1:length(nodeList)){
  eval(parse(text = paste("if(length(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]) > 0){floodedNodesMaster$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "[which(floodedNodesMaster$Node == \"", nodeList[k], "\")] <- nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodDuration_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]}", sep = "")))
  eval(parse(text = paste("if(length(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]) > 0){floodedNodesMaster$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "[which(floodedNodesMaster$Node == \"", nodeList[k], "\")] <- nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$FloodVolume_", modelTypes[i], "_", scenarioNames[j], "[which(nodeFlooding_", modelTypes[i], "_", scenarioNames[j], "$Node == \"", nodeList[k], "\")]}", sep = "")))
  
}


##These get rid of some extraneous columns that exist in the base file
floodedNodesMaster <- floodedNodesMaster[,-3]
floodedNodesMaster <- floodedNodesMaster[,-2]
floodedNodesMaster <- data.frame(floodedNodesMaster)

##Writing the flooded nodes table generated in the above code to a csv file that was later joined to the shapefile in ArcGIS
write.csv(floodedNodesMaster, file = "floodedNodesMaster_R.csv")


#Creating plots of node flooding trends based on scenario and model type, Also Generating p values and r-squared values for all scenarios and model types to determine which will get trendlines
nodeFloodingInformation <- data.frame(read.csv("NodeFloodingInformation.csv"))
varNames <- colnames(nodeFloodingInformation)[c(4,5,7,8,10)]
varNamesForColumns <- c(paste(varNames, "_pValue", sep = ""), paste(varNames, "_adjRSquared", sep = ""))
featureType <- unique(nodeFloodingInformation$Feature)
pValueFrameNodeFlooding <- data.frame(matrix(nrow = 3, ncol = 11, 0))
colnames(pValueFrameNodeFlooding) <- c("Feature", varNamesForColumns)
pValueFrameNodeFlooding$Feature <- featureType


##Calculating p-values and r-squareds. Running for all three model builds and for all variables
for(i in 1:length(featureType)){
  for(j in 1:length(varNames)){
    eval(parse(text = paste("yVariables <- nodeFloodingInformation[which(nodeFloodingInformation$Feature == \"", featureType[i], "\"),]", sep = "")))
    eval(parse(text = paste("lm_", featureType[i], "_", varNames[j], " <- lm(", varNames[j], " ~ WetlandLoss, yVariables)", sep = "")))
    eval(parse(text = paste("pValueFrameNodeFlooding[i, 1 + j] <- glance(lm_", featureType[i], "_", varNames[j], ")$p.value", sep = "")))
    eval(parse(text = paste("pValueFrameNodeFlooding[i, 6 + j] <- summary(lm_", featureType[i], "_", varNames[j], ")$adj.r.squared", sep = "")))
  }
}

write.csv(pValueFrameNodeFlooding, file = "nodeFloodingPValues.csv")


##Plotting node flooding trends

g_FloodNodes <- ggplot(data = nodeFloodingInformation, aes(x = WetlandLoss, y = NumberFloodedNodesAllSystem)) +
  geom_point(aes(shape = Feature), size = 3)+
  geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black")+
  theme_classic()+
  theme(legend.position = "none")



g_FloodNodes_1 <- g_FloodNodes + labs(title = "Number of nodes flooded",
                                      y = "Number of flooded nodes",
                                      x = "Lost wetland area", color = "") +
  scale_linetype_manual(values = c("dashed", "solid", "blank")) +
  scale_shape_manual(values = c(1, 17, 0)) + 
  scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))



g_WetlandNodeFlooding <- ggplot(data = nodeFloodingInformation, aes(x = WetlandLoss, y = FloodVolumeInWetlandNodes_1000m3)) +
  geom_point(aes(shape = Feature), size = 3)+
  geom_smooth(aes(linetype = Feature), method = "lm", se = FALSE, color = "black")+
  theme_classic()+
  theme(legend.position = "none")

g_WetlandNodeFlooding_1 <- g_WetlandNodeFlooding + labs(title = "Flooding from wetland nodes",
                                                        y = expression(paste("Flood volume (1000 ", m^3, ")", sep = "")),
                                                        x = "Lost wetland area", color = "") +
  scale_linetype_manual(values = c("dashed", "solid", "dotted")) +
  scale_shape_manual(values = c(1, 17, 0)) + 
  scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))

g_PercentWetlandFlood <- ggplot(data = nodeFloodingInformation, aes(x = WetlandLoss, y = PercentFloodVolumeByWetlandNodes)) +
  geom_point(aes(shape = Feature), size = 3)+
  geom_smooth(aes(linetype = Feature), method = "lm", se = FALSE, color = "black")+
  theme_classic()+
  theme(legend.position = "none")

g_PercentWetlandFlood_1 <- g_PercentWetlandFlood + labs(title = "Wetland contribution to flooding",
                                                        y = "Proportion of flood volume from wetlands",
                                                        x = "Lost wetland area", color = "") +
  scale_linetype_manual(values = c("dashed", "solid", "blank")) +
  scale_shape_manual(values = c(1, 17, 0)) + 
  scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))

g_WetlandNodesSpentFlooding <- ggplot(data = nodeFloodingInformation, aes(x = WetlandLoss, y = FloodTimeInWetlandNodes_hours)) +
  geom_point(aes(shape = Feature), size = 3)+
  geom_smooth(aes(linetype = Feature), method = "lm", se = FALSE, color = "black")+
  theme_classic()+
  theme(legend.position = "none")

g_WetlandNodesSpentFlooding_1 <- g_WetlandNodesSpentFlooding + labs(title = "Time wetlands spent flooding",
                                                                    y = "Flood duration (hours)",
                                                                    x = "Lost wetland area", color = "") +
  scale_linetype_manual(values = c("dashed", "solid", "blank")) +
  scale_shape_manual(values = c(1, 17, 0)) + 
  scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))


g_WetlandFloodTimeContribution <- ggplot(data = nodeFloodingInformation, aes(x = WetlandLoss, y = PercentWetlandFloodTimeContribution)) +
  geom_point(aes(shape = Feature), size = 3)+
  geom_smooth(aes(linetype = Feature), method = "lm", se = FALSE, color = "black")+
  theme_classic()+
  theme(legend.position = "none")

g_WetlandFloodTimeContribution_1 <- g_WetlandFloodTimeContribution + labs(title = "Proportion of time wetlands contribute to flooding",
                                                                          y = "Proportion of flood time from wetlands",
                                                                          x = "Lost wetland area", color = "") +
  scale_linetype_manual(values = c("dashed", "solid", "blank")) +
  scale_shape_manual(values = c(1, 17, 0)) + 
  scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))



pdf(file = "NodeFloodingLocations_02022021.pdf")
ggarrange(g_FloodNodes_1, g_WetlandNodeFlooding_1, 
          g_WetlandNodesSpentFlooding_1, g_WetlandFloodTimeContribution_1, 
          g_PercentWetlandFlood_1,
          ncol = 2, nrow =2)
dev.off()


#Plotting downstream flooding of wetland 43 (also referred to as Wetland 11) to produce the figure that demonstrates how flooding is tabulated.
floodWetland11 <- data.frame(read.csv("Wetland11FloodingDepth.csv"))

g_FloodWetland11 <- ggplot(floodWetland11, aes (x = Hour, y = Flooding))+
  geom_line(linetype = "solid", color = "black") +
  theme_classic() +
  ylab("Rate of flooding (m /s)") +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

g_DepthWetland11 <- ggplot(floodWetland11, aes (x = Hour, y = Depth))+
  geom_line(linetype = "solid", color = "black") +
  theme_classic() +
  ylab("Depth (m)") +
  theme(legend.position = "none",) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

gHyetograph <- ggplot(precipitationForHyetograph, aes(Hour, Precipitation)) +
  geom_bar(stat = "identity", fill = "black") +
  theme_classic() +
  ylab("Precip") +
  labs(title = "") +
  scale_y_reverse(labels = scales::number_format(accuracy = 0.1)) +
  #scale_x_continuous(labels = scales::number_format(accuracy = 0.1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank())


pdf(file = "PlotNodeFloodingDepthWetland11.pdf")
grid.arrange(gHyetograph, g_FloodWetland11, g_DepthWetland11, ncol = 1, heights = c(1,3,3))
dev.off()



#Creating tables for each of the sectors, to show how wetland loss changes flood volume and location and stress on network
resultsActual <- data.frame(read.csv("D:/Valdivia/DissertationChapter3/ResultsActual.csv"), stringsAsFactors = FALSE)
resultsBarriosBajos <- data.frame(read.csv("D:/Valdivia/DissertationChapter3/ResultsBarriosBajos.csv"), stringsAsFactors = FALSE)
resultsTeja <- data.frame(read.csv("D:/Valdivia/DissertationChapter3/ResultsTeja.csv"), stringsAsFactors = FALSE)
resultsNorte <- data.frame(read.csv("D:/Valdivia/DissertationChapter3/ResultsNorte.csv"), stringsAsFactors = FALSE)

##Expanding results Norte data frame so that it can be combined with the 
##others more easily in the for loop
resultsNorte <- rbind(resultsNorte, resultsNorte, resultsNorte)

##Creating data frame to sum all the flooding data across all sectors
resultsSum <- as.data.frame(matrix(nrow = 18, ncol = 10, 0))
colnames(resultsSum) <- colnames(resultsActual)

resultsSum[,1:3] <- resultsActual[,1:3]
resultsSum[,10] <- resultsActual[,10]
sectorNames <- c("Actual", "BarriosBajos", "Norte", "Teja")

##Filling the resultsSum data frame using the data included in the other results
##files

for(i in 1:nrow(resultsActual)){
  for(j in 4:(ncol(resultsActual)-1)){
    for(k in 1:length(sectorNames)){
      eval(parse(text = paste("resultsSum[i,j] <- resultsSum[i,j] + results", sectorNames[k], "[i,j]", sep = "")))
    }    
  }
}

resultsActual <- data.frame(read.csv("ResultsActual_11102020.csv"))
resultsNorte <- data.frame(read.csv("ResultsNorte_11102020.csv"))
resultsTeja <- data.frame(read.csv("ResultsTeja_11102020.csv"))
resultsBarriosBajos <- data.frame(read.csv("ResultsBarriosBajos_11102020.csv"))
resultsSum <- data.frame(read.csv("ResultsAllSectors_11102020.csv"))


##Note, csvs originally had the incorrect wetland area lost attached to them, but that has been rectified. This step is unnecessary now, but I'm leaving it in in case later code uses it
lostWetlandArea <- c(0.3731, 0.2388, 0.1827, 0.1326, 0.0972, 0.0000)
resultsSum$WetlandLoss <- lostWetlandArea
resultsSum <- resultsSum[,2:11]




#Collecting significance for all the system-wide wetland loss variables into data frame
##Building variables and data frame
varNames <- colnames(resultsSum)[4:9]
varNamesForColumns <- c(paste(varNames, "_pValue", sep = ""), paste(varNames, "_adjRSquared", sep = ""))
featureType <- unique(resultsSum$Feature)
pValueFrame <- data.frame(matrix(nrow = 3, ncol = 13, 0))
colnames(pValueFrame) <- c("Feature", varNamesForColumns)
pValueFrame$Feature <- featureType
i = 1
j = 1

##Running for all three model builds and for all variables
for(i in 1:length(featureType)){
  for(j in 1:length(varNames)){
    eval(parse(text = paste("yVariables <- resultsSum[which(resultsSum$Feature == \"", featureType[i], "\"),]", sep = "")))
    eval(parse(text = paste("lm_", featureType[i], "_", varNames[j], " <- lm(", varNames[j], " ~ WetlandLoss, yVariables)", sep = "")))
    eval(parse(text = paste("pValueFrame[i, 1 + j] <- glance(lm_", featureType[i], "_", varNames[j], ")$p.value", sep = "")))
    eval(parse(text = paste("pValueFrame[i, 7 + j] <- summary(lm_", featureType[i], "_", varNames[j], ")$adj.r.squared", sep = "")))
  }
}

write.csv(pValueFrame, file = "wetlandLossPValues.csv")



#Plotting system-wide wetland loss effects for all model types (wet, wshed, wet + wshed)
yNames <- c("expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\"))", "expression(paste(\"Hours nodes flooded\", sep = \"\"))", "expression(paste(\"Hours nodes surcharged\", sep = \"\"))", "expression(paste(\"Max stored volume (1000 \", m^3, \")\", sep = \"\"))", "expression(paste(\"Hours capacity limited\", sep = \"\"))", "expression(paste(\"Hours above normal\", sep = \"\"))")
titleNames <- c("Node flooding", "Hours nodes flooded", "Hours nodes surcharged", "Wetland storage volume", "Hours conduit capacity limited", "Hours flow above normal")
varNames <- colnames(resultsSum)[4:9]


eval(parse(text = paste("g_", varNames, " <- ggplot(data = resultsSum, aes(x = WetlandLoss, y = ", varNames, ", group = Feature)) +
                        labs(y = ", yNames, ") +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")",
                        sep = "")))


eval(parse(text = paste("g_", varNames, "_1 <- g_", varNames, " + labs(title = \"", titleNames, "\",
                        y = ", yNames, ",
                        x = \"Lost wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"solid\", \"dotted\")) + 
                        scale_shape_manual(values = c(1, 17, 0)) +
                        scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))", sep = "")))

#Getting rid of trendline for relevant plots
g_NodeFloodVolume_Wetlands <- ggplot(data = resultsSum, aes(x = WetlandLoss, y = NodeFloodVolume_Wetlands, group = Feature)) +
                        labs(y = expression(paste("Flood volume (1000 ", m^3, ")", sep = ""))) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black")+
                        theme_classic()+
                        theme(legend.position = "none")


g_NodeFloodVolume_Wetlands_1 <- g_NodeFloodVolume_Wetlands + labs(title = "Node Flooding",
                        y = expression(paste("Flood volume (1000 ", m^3, ")", sep = "")),
                        x = "Lost wetland area", color = "") +
                        scale_linetype_manual(values = c("dashed", "solid", "dotted")) + 
                        scale_shape_manual(values = c(1, 17, 0)) +
                        scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))


g_HoursNodesFlooded_Wetlands <- ggplot(data = resultsSum, aes(x = WetlandLoss, y = HoursNodesFlooded_Wetlands, group = Feature)) +
                        labs(y = "") +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black")+
                        theme_classic()+
                        theme(legend.position = "none")


g_HoursNodesFlooded_Wetlands_1 <- g_HoursNodesFlooded_Wetlands + labs(title = "Hours nodes flooded",
                        y = "Duration of flooding (hours)",
                        x = "Lost wetland area", color = "") +
                        scale_linetype_manual(values = c("dashed", "solid", "blank")) + 
                        scale_shape_manual(values = c(1, 17, 0)) +
                        scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))


g_MaxStoredVolume_Wetlands <- ggplot(data = resultsSum, aes(x = WetlandLoss, y = MaxStoredVolume_Wetlands, group = Feature)) +
                        labs(y = "") +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black")+
                        theme_classic()+
                        theme(legend.position = "none")


g_MaxStoredVolume_Wetlands_1 <- g_MaxStoredVolume_Wetlands + labs(title = "Max stored volume",
                        y = expression(paste("Stored volume (1000 ", m^3, ")", sep = "")),
                        x = "Lost wetland area", color = "") +
                        scale_linetype_manual(values = c("dashed", "solid", "blank")) + 
                        scale_shape_manual(values = c(1, 17, 0)) +
                        scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))


g_HoursAboveNormalFlow_Wetlands <- ggplot(data = resultsSum, aes(x = WetlandLoss, y = HoursAboveNormalFlow_Wetlands, group = Feature)) +
                        labs(y = "Time above normal flow (hours)") +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black")+
                        theme_classic()+
                        theme(legend.position = "none")


g_HoursAboveNormalFlow_Wetlands_1 <- g_HoursAboveNormalFlow_Wetlands + labs(title = "Hours above normal flow",
                        y = "Hours above normal flow",
                        x = "Lost wetland area", color = "") +
                        scale_linetype_manual(values = c("dashed", "solid", "blank")) + 
                        scale_shape_manual(values = c(1, 17, 0)) +
                        scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))


g_HoursConduitCapacityLimited_Wetlands <- ggplot(data = resultsSum, aes(x = WetlandLoss, y = HoursConduitCapacityLimited_Wetlands, group = Feature)) +
                        labs(y = expression(paste("Flood volume (1000", m^3, ")", sep = ""))) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black")+
                        theme_classic()+
                        theme(legend.position = "none")


g_HoursConduitCapacityLimited_Wetlands_1 <- g_HoursConduitCapacityLimited_Wetlands + labs(title = "Hours conduit capacity limited",
                        y = "Time conduit capacity limited (hours)",
                        x = "Lost wetland area", color = "") +
                        scale_linetype_manual(values = c("dashed", "solid", "blank")) + 
                        scale_shape_manual(values = c(1, 17, 0)) +
                        scale_x_continuous(limits = c(0, 0.42), breaks = seq(0, 0.50, 0.10), labels = scales::percent(seq(0, 50, 10), scale = 1))





#All sectors
pdf(file = "D:/Valdivia/DissertationChapter3/RCode/LostWetlandAreaGraphics_AllSectors_01252021.pdf")
ggarrange(g_NodeFloodVolume_Wetlands_1, g_HoursNodesFlooded_Wetlands_1,
          ncol = 2, nrow =2)

ggarrange(g_HoursNodesSurcharged_Wetlands_1, g_MaxStoredVolume_Wetlands_1,
          ncol = 2, nrow =2)

ggarrange(g_HoursConduitCapacityLimited_Wetlands_1, g_HoursAboveNormalFlow_Wetlands_1,
          ncol = 2, nrow =2)

dev.off()



#Writing results tables
write.csv(resultsActual, "ResultsActual_11102020.csv")
write.csv(resultsNortePlots, "ResultsNorte_11102020.csv")
write.csv(resultsTeja, "ResultsTeja_11102020.csv")
write.csv(resultsBarriosBajos, "ResultsBarriosBajos_11102020.csv")
write.csv(resultsSum, "ResultsAllSectors_11102020.csv")


#Plots for relating wetland 43 area loss and downstream node flooding
resultsWetland43_downDrain <- read.csv("D:/Valdivia/DissertationChapter3/ResultsWetland43_downDrain.csv", stringsAsFactors = FALSE)
resultsWetland43 <- read.csv("D:/Valdivia/Dissertationchapter3/ResultsWetland43.csv")
resultsWetland43_downDrain2 <- read.csv("D:/Valdivia/DissertationChapter3/ResultsWetland43_downDrain2.csv", stringsAsFactors = FALSE)

#Relating the results from the models where ponding was allowed to the models where ponding was not allowed. Providing aggregate results instead of node-by-node results for simplification.
scenarios <- c("Business-as-usual", "Eco-wetland", "Friendly", "Inclusive", "Resilient to flooding", "Present Day")

resultsWetland43$Feature[which(resultsWetland43$Feature == "Wet_Wshed")] <- "WetlandWshed"
resultsWetland43_downDrain$Feature[which(resultsWetland43_downDrain$Feature == "Wet_Wshed")] <- "WetlandWshed"
resultsWetland43_downDrain2$Feature[which(resultsWetland43_downDrain2$Feature == "Wet_Wshed")] <- "WetlandWshed"

resultsWetland43$Feature[which(resultsWetland43$Feature == "Wet")] <- "Wetland"
resultsWetland43_downDrain$Feature[which(resultsWetland43_downDrain$Feature == "Wet")] <- "Wetland"
resultsWetland43_downDrain2$Feature[which(resultsWetland43_downDrain2$Feature == "Wet")] <- "Wetland"

floodVolumes_Master$Scenario[which(floodVolumes_Master$Scenario == "Present day")] <- "Present Day"

resultsWetland43_difference <- data.frame(matrix(nrow = 0, ncol = 6, NA))
colnames(resultsWetland43_difference) <- c("Scenario", "Feature", "WetlandLoss", "NodeFlooding_Ponding", "NodeFlooding_noPonding", "NodeFlooding_Difference" )
i = 1
j = 1

for(i in 1:length(scenarios)){
  for(j in 1:length(modelType)){
    eval(parse(text = paste("flooding_noPonding <- sum(resultsWetland43$NodeFloodVolume_Wetlands[which(resultsWetland43$Scenario == \"", scenarios[i], "\" & resultsWetland43$Feature == \"", modelType[j], "\")], resultsWetland43_downDrain$NodeFloodVolume_Wetlands[which(resultsWetland43_downDrain$Scenario == \"", scenarios[i], "\" & resultsWetland43_downDrain$Feature == \"", modelType[j], "\")], resultsWetland43_downDrain2$NodeFloodVolume_Wetlands[which(resultsWetland43_downDrain2$Scenario == \"", scenarios[i], "\" & resultsWetland43_downDrain2$Feature == \"", modelType[j], "\")])", sep = "")))
    eval(parse(text = paste("flooding_ponding <- sum(floodVolumes_Master$NodeFlooding[which(floodVolumes_Master$Scenario == \"", scenarios[i], "\" & floodVolumes_Master$Feature == \"", modelType[j], "\")])", sep = "")))
    flooding_difference <- flooding_ponding - flooding_noPonding
    floodRow <- data.frame(matrix(nrow = 1, ncol = 6, NA))
    colnames(floodRow) <- c("Scenario", "Feature", "WetlandLoss", "NodeFlooding_Ponding", "NodeFlooding_noPonding", "NodeFlooding_Difference" )
    eval(parse(text = paste("floodRow[1,1] <- \"", scenarios[i], "\"", sep = "")))
    eval(parse(text = paste("floodRow[1,2] <- \"", modelType[j], "\"", sep = "")))
    if(scenarios[i] == "Business-as-usual"){floodRow[1,3] <- 0.97} else if(scenarios[i] == "Eco-wetland"){floodRow[1,3] <- 0.76} else if(scenarios[i] == "Friendly"){floodRow[1,3] <- 0.30} else if(scenarios[i] == "Inclusive"){floodRow[1,3] <- 0.35} else if(scenarios[i] == "Resilient to flooding"){floodRow[1,3] <- 0.66} else if(scenarios[i] == "Present Day"){floodRow[1,3] <- 0.00}
    floodRow[1,4] <- flooding_ponding
    floodRow[1,5] <- flooding_noPonding
    floodRow[1,6] <- flooding_difference
    resultsWetland43_difference <- rbind(resultsWetland43_difference, floodRow)
  }
}

#write.csv(resultsWetland43_difference, "resultsWetland43_difference.csv")

#Seeing if there are significant trends
wetland43_difference_pValues <- data.frame(matrix(ncol = 3, nrow = 0, NA))
colnames(wetland43_difference_pValues) <- c("modelType", "pValue", "rSquared")

j = 1
j = 2
for(j in 1:length(modelType)){
  eval(parse(text = paste(modelType[j], "_pValues <- data.frame(matrix(ncol = 3, nrow = 1, NA))", sep = "")))
  eval(parse(text = paste("colnames(", modelType[j], "_pValues) <- c(\"modelType\", \"pValue\", \"rSquared\")", sep = "")))
  eval(parse(text = paste(modelType[j], "_pValues$modelType <- \"", modelType[j], "\"", sep = "")))
  eval(parse(text = paste("floodVolumes_", modelType[j], " <- resultsWetland43_difference[which(resultsWetland43_difference$Feature == \"", modelType[j], "\"),]", sep = "")))
  eval(parse(text = paste("lm_floodVolumes_", modelType[j], " <- lm(NodeFlooding_Difference ~ WetlandLoss, floodVolumes_", modelType[j], ")", sep = "")))
  eval(parse(text = paste(modelType[j], "_pValues$pValue <- glance(lm_floodVolumes_", modelType[j], ")$p.value", sep = "")))
  eval(parse(text = paste(modelType[j], "_pValues$rSquared <- summary(lm_floodVolumes_", modelType[j], ")$adj.r.squared", sep = "")))
}
eval(parse(text = paste("wetland43_difference_pValues <- rbind(wetland43_difference_pValues, ", modelType, "_pValues)", sep = "")))

#write.csv(wetland43_difference_pValues, "wetland43_difference_pValues.csv")


g_wet43_difference <- ggplot(data = resultsWetland43_difference, aes(x = WetlandLoss, y = NodeFlooding_Difference, group = Feature)) +
  geom_point(aes(shape = Feature), size = 3) +
  geom_smooth(method = "lm", aes(linetype = Feature), se = FALSE, color = "black") +
  theme_classic() +
  theme(legend.position = "none")

g_wet43_difference_1 <- g_wet43_difference + labs(title = "Model comparison", y = expression(paste("Flood volume difference (1000 ", m^3, ")", sep = "")), 
                                                  x = "Lost wetland area", color = "") +
  scale_linetype_manual(values = c("solid", "dashed", "blank")) +
  scale_shape_manual(values = c(1, 17, 0)) +
  scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1))
  

pdf(file = "wet43FloodVolumeDifference_04282021.pdf")
ggarrange(g_wet43_difference_1,
          ncol = 1, nrow = 1)
dev.off()


pdf(file = "D:/Valdivia/DissertationChapter3/RCode/LostWetlandAreaGraphics_BarriosBajos.pdf")
ggarrange(g_NodeFloodVolume_Wetlands_BB_1, g_HoursNodesFlooded_Wetlands_BB_1,
          labels = c("A", "B"),
          ncol = 1, nrow =2)

ggarrange(g_HoursNodesSurcharged_Wetlands_BB_1, g_MaxStoredVolume_Wetlands_BB_1,
          labels = c("C", "D"),
          ncol = 1, nrow =2)

ggarrange(g_HoursConduitCapacityLimited_Wetlands_BB_1, g_HoursAboveNormalFlow_Wetlands_BB_1,
          labels = c("E", "F"),
          ncol = 1, nrow =2)

dev.off()


#titleNames <- c("Node flood volume vs wetland loss", "Hours nodes flooded vs wetland loss", "Hours nodes surcharged vs wetland loss", "Max stored volume vs wetland loss", "Hours conduit capacity limited vs wetland loss", "Hours flow above normal vs wetland loss")
#yAxesLabels <- c("Flood volume (1000 m^3)", "Hours nodes flooded", "Hours nodes surcharged", "Max stored volume (1000 m3)", "Hours conduit capacity limited", "Hours flow above normal") 
eval(parse(text = paste("g_Wet43 <- ggplot(data = resultsWetland43, aes(x = WetlandLoss, y = NodeFloodVolume_Wetlands, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))

g_Wet43

eval(parse(text = paste("g_Wet43_1 <- g_Wet43 + labs(title = \"Updrain wetland\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"solid\", \"blank\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1))",
                        sep = "")))

g_Wet43_1
eval(parse(text = paste("g_Wet43_downDrain <- ggplot(data = resultsWetland43_downDrain, aes(x = WetlandLoss, y = NodeFloodVolume_Wetlands, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))


eval(parse(text = paste("g_Wet43_downDrain_1 <- g_Wet43_downDrain + labs(title = \"Downdrain node 1\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost updrain wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"solid\", \"dotted\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1))",
                        sep = "")))
g_Wet43_downDrain_1


eval(parse(text = paste("g_Wet43_downDrain2 <- ggplot(data = resultsWetland43_downDrain2, aes(x = WetlandLoss, y = NodeFloodVolume_Wetlands, group = Feature)) +
                        geom_point(aes(shape = Feature), size = 3)+
                        geom_smooth(method = \"lm\", aes(linetype = Feature), se = FALSE, color = \"black\")+
                        theme_classic()+
                        theme(legend.position = \"none\")", sep = "")))


eval(parse(text = paste("g_Wet43_downDrain2_1 <- g_Wet43_downDrain2 + labs(title = \"Downdrain node 2\",
                        y = expression(paste(\"Flood volume (1000 \", m^3, \")\", sep = \"\")),
                        x = \"Lost updrain wetland area\", color = \"\") +
                        scale_linetype_manual(values = c(\"dashed\", \"blank\", \"blank\")) +
                        scale_shape_manual(values = c(1, 17, 0)) + 
                        scale_x_continuous(limits = c(0, 1.02), breaks = seq(0, 1, 0.20), labels = scales::percent(seq(0, 100, 20), scale = 1))",
                        sep = "")))

g_Wet43_downDrain2_1


pdf(file = "Wetland43_LinearPlot_01202021.pdf")
ggarrange(g_Wet43_1, g_Wet43_1, g_Wet43_downDrain_1, g_Wet43_downDrain2_1,
          ncol = 2, nrow =2)
dev.off()
